/* Functions for Marquardt-Levenberg minimization routine.
   All functions are of the form y = F(x, p, dydp, npar).

   By:   S.C. Molitor (smolitor@med.unc.edu)
   Date: February 29, 2000 */


// inclusions

#include <math.h>
#include "mex.h"


// EXPONENTIAL - a sum of N exponential functions:
//
//    F(X) = A0 + A1*exp(-t/T1) + A2*exp(-t/T2) + ... + AN*exp(-t/TN)
//
// The parameter vector P must have 2N + 1 parameters for N exponentials:
//
//                 P = [A0 A1 T1 A2 T2 ... AN TN]

double exponential(double	x,		// x value, returns f(p, x)
				   double	p[],	// parameter values in p
				   double	dfdp[],	// partial derivatives df/dp
				   int		npar) {	// number of parameters in p

	// declare variables
	
	bool	errflag;	// error flag
	double	y, yi;		// evaluated f(p, x)
	int		nexp;		// number of exponentials added
	int		i;			// loop counter

	// calculate number of exponentials added together
	// check number of parameters
	// check for positive time constants

	errflag = false;
	nexp = (npar - 1) / 2;
	if (npar == (2 * nexp + 1)) {
		for (i = 1; i <= nexp; i++) {
			if (p[2 * i] <= 0.0) {
				errflag = true;
				break;
			}
		}
	}
	else {
		errflag = true;
	}

	// if error, return zero for f(p, x) & df/dp

	if (errflag) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// calculate df/dp
	// sum individual exponentials

	y = p[0];
	dfdp[0] = 1.0;
	for (i = 1; i <= nexp; i++) {
		dfdp[2 * i - 1] = exp(-x / p[2 * i]);
		yi = p[2 * i - 1] * dfdp[2 * i - 1];
		dfdp[2 * i] = yi * x / (p[2 * i] * p[2 * i]);
		y += yi;
	}
	return(y);
}


// EXPPOWER - an exponential function raised to the Nth power:
//
//           F(X) = Imax*[1 + (A0 - 1)*exp(-t/T)]^N
//
// where 0 <= A0 <= 1.  The parameter vector P must have 4
// parameters:
//                     P = [Imax A0 T N]

double exppower(double	x,		// x value, returns f(p, x)
			    double	p[],	// parameter values in p
			    double	dfdp[],	// partial derivatives df/dp
				int		npar) {	// number of parameters in p

	// declare variables
	
	bool	errflag;	// error flag
	double	amp;		// amplitude (A0 - 1)
	double	expfac;		// exponential factor
	double	exparg;		// exponential argument
	double	expargn_1;	// exponential argument ^ n - 1
	double	y;			// evaluated f(p, x)
	int		i;			// loop counter

	// check number of parameters
	// check for 0 <= A0 <= 1
	// check for positive time constant
	// check for N >= 1

	if (npar != 4) {
		errflag = true;
	}
	else if ((p[1] < 0.0) || (p[1] > 1.0)) {
		errflag = true;
	}
	else if (p[2] <= 0.0) {
		errflag = true;
	}
	else if (p[3] < 1) {
		errflag = true;
	}
	else {
		errflag = false;
	}

	// if error, return zero for f(p, x) & df/dp

	if (errflag) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// calculate values for f(x) & df/dp
	// exponential arg should be > 0 when 0 <= A0 & A1 <= 1
	// return zero for f(p, x) & df/dp if not the case

	amp = p[1] - 1.0;
	expfac = exp(-x / p[2]);
	exparg = 1.0 + amp * expfac;
	if (exparg > 0.0) {
		expargn_1 = pow(exparg, p[3] - 1.0);
	}
	else {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// calculate f(x)
	// calculate df/dp

	dfdp[0] = exparg * expargn_1;
	y = p[0] * dfdp[0];
	dfdp[1] = p[0] * expargn_1 * p[3] * expfac;
	dfdp[2] = p[0] * expargn_1 * p[3] * amp * x * expfac / (p[2] * p[2]);
	dfdp[3] = y * log(exparg);
	return(y);
}


// GAUSSIAN - a sum of N scaled Gaussian PDF functions:
//
//        F(X) = A0 + A1*G1(X) + A2*G2(X) + ... + AN*GN(X)
//
// where each Gaussian Gi(X) has a mean Mi and a standard deviation Si:
//
//     Gi(X) = exp(-(X - Mi)^2 / (2*Si^2)) / (SQRT(2*PI) * Si)
//
// The parameter vector P must have 3N + 1 parameters for N Gaussians:
//
//            P = [A0 A1 M1 S1 A2 M2 S2 ... AN MN SN]

double gaussian(double	x,		// x value, returns f(p, x)
				double	p[],	// parameter values in p
				double	dfdp[],	// partial derivatives df/dp
				int		npar) {	// number of parameters in p

	// declare variables
	
	bool	errflag;	// error flag
	double	ampfac;		// amplitude scaling
	double	normx;		// x normalized by mu & omega
	double	expfac;		// exponential factor
	double	y, yi;		// evaluated f(p, x)
	int		npdf;		// number of gaussian PDFs added
	int		i;			// loop counter

	// calculate number of Gaussian PDFs added together
	// check number of parameters
	// check for positive std. deviations

	errflag = false;
	npdf = (npar - 1) / 3;
	if (npar == (3 * npdf + 1)) {
		for (i = 1; i <= npdf; i++) {
			if (p[3 * i] <= 0.0) {
				errflag = true;
				break;
			}
		}
	}
	else {
		errflag = true;
	}

	// if error, return zero for f(p, x) & df/dp

	if (errflag) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// calculate df/dp
	// sum individual Gaussian PDFs

	y = p[0];
	dfdp[0] = 1.0;
	for (i = 1; i <= npdf; i++) {
		ampfac = 1.0 / (sqrt(8.0 * atan2(1, 1)) * p[3 * i]);
		normx = (x - p[3 * i - 1]) / p[3 * i];
		expfac = exp(-0.5 * normx * normx);
		dfdp[3 * i - 2] = ampfac * expfac;
		yi = p[3 * i - 2] * dfdp[3 * i - 2];
		dfdp[3 * i - 1] = yi * normx / p[3 * i];
		dfdp[3 * i] = yi * (normx * normx - 1.0) / p[3 * i];
		y += yi;
	}
	return(y);
}


// BOLTZMANN - multistate Boltzmann function:
//
//      F(X) = A0 + A1/(1 + EN + EN*EN-1 + EN*EN-1*EN-2 + ...)
//
// where each Boltzmann Ei(X) has a slope Ki and a half-activation Xi:
//
//                   Ei(X) = exp(-Ki*(X - Xi))
//
// The parameter vector P must have 2N + 2 elements for N states:
//
//               P = [A0 A1 K1 X1 K2 X2 ... KN XN]

double boltzmann(double	x,		// x value, returns f(p, x)
				 double	p[],	// parameter values in p
				 double	dfdp[],	// partial derivatives df/dp
				 int	npar) {	// number of parameters in p

	// declare variables
	
	double	product;	// product of Boltzmanns
	double	expfac;		// exponential calculation
	double	*dfdpfac;	// storage for df/dp scaling
	double	y;			// evaluated f(p, x)
	int		nboltz;		// number of Boltzmanns added
	int		i, j;		// loop counters

	// calculate number of Boltzmanns added together
	// check number of parameters
	// if error, return zero for f(p, x) & df/dp

	nboltz = (npar - 2) / 2;
	if (npar != (2 * nboltz + 2)) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// allocate space for df/dp calculations

	dfdpfac = (double *) mxCalloc(nboltz, sizeof(double));

	// calculate individual Boltzmann values
	// calculate sums & product of Boltzmanns

	for (product = 1.0, i = nboltz; i >= 1; i--) {
		expfac = exp(-p[2 * i] * (x - p[2 * i + 1]));
		product *= expfac;
		for (j = i - 1; j < nboltz; j++) {
			dfdpfac[j] += product;
		}
	}

	// complete calculations for f(p, x) and df/dp
	// free allocated space
	// return f(p, x)

	dfdp[0] = 1.0;
	dfdp[1] = 1.0 / (1.0 + dfdpfac[nboltz - 1]);
	y = p[0] + p[1] * dfdp[1];
	expfac = p[1] * dfdp[1] * dfdp[1];
	for (i = 1; i <= nboltz; i++) {
		dfdp[2 * i] = expfac * dfdpfac[i - 1] * (x - p[2 * i + 1]);
		dfdp[2 * i + 1] = expfac * dfdpfac[i - 1] * (-p[2 * i]);
	}
	mxFree((void *) dfdpfac);
	return(y);
}


// BOLTZSUM - sum of N individual Boltzmann functions:
//
//     F(X) = A0 + A1/(1 + E1) + A2/(1 + E2) + ... + AN/(1 + EN)
//
// where each Boltzmann Ei(X) has a slope Ki and a half-activation Xi:
//
//                   Ei(X) = exp(-Ki*(X - Xi))
//
// The parameter vector P must have 3N + 1 elements for N individual
// Boltzmanns:
//
//            P = [A0 A1 K1 X1 A2 K2 X2 ... AN KN XN]

double boltzsum(double	x,		// x value, returns f(p, x)
			    double	p[],	// parameter values in p
			    double	dfdp[],	// partial derivatives df/dp
				int		npar) {	// number of parameters in p

	// declare variables
	
	double	expfac;		// exponential calculation (Ei)
	double	expprod;	// exponential calculations (1/(1 + Ei) & Ei/(1 + Ei)^2)
	double	y;			// evaluated f(p, x)
	int		nboltz;		// number of Boltzmanns added
	int		i;			// loop counter

	// calculate number of Boltzmanns added together
	// check number of parameters
	// if error, return zero for f(p, x) & df/dp

	nboltz = (npar - 1) / 3;
	if (npar != (3 * nboltz + 1)) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// initialize f(x) & df/dp

	y = p[0];
	dfdp[0] = 1.0;

	// loop through Boltzmann sums
	// add sums to f(x)
	// calculate relevant df/dp at each sum

	for (i = 1; i <= nboltz; i++) {
		expfac = exp(-p[3 * i - 1] * (x - p[3 * i]));
		expprod = 1.0 / (1.0 + expfac);
		dfdp[3 * i - 2] = expprod;
		y += p[3 * i - 2] * expprod;
		expprod *= expprod * expfac;
		dfdp[3 * i - 1] = p[3 * i - 2] * (x - p[3 * i]) * expprod;
		dfdp[3 * i] = -p[3 * i - 2] * p[3 * i - 1] * expprod;
	}
	return(y);
}


// IVBOLTZSUM - current form for sum of N Boltzmann functions:
//
// F(X) = A0 + (X - Vr)*[A1/(1 + E1) + A2/(1 + E2) + ... + AN/(1 + EN)]
// 
// where each Boltzmann Ei(X) has a slope Ki and a half-activation Xi:
//
//                   Ei(X) = exp(-Ki*(X - Xi))
//
// The parameter vector P must have 3N + 2 elements for N individual
// Boltzmanns:
//
//           P = [A0 Vr A1 K1 X1 A2 K2 X2 ... AN KN XN]
//

double ivboltzsum(double	x,		// x value, returns f(p, x)
			      double	p[],	// parameter values in p
			      double	dfdp[],	// partial derivatives df/dp
				  int		npar) {	// number of parameters in p

	// declare variables
	
	double	drvforce;	// driving force (x - Vr)
	double	expfac;		// exponential calculation (Ei)
	double	expprod;	// exponential calculations (1/(1 + Ei) & Ei/(1 + Ei)^2)
	double	y;			// evaluated f(p, x)
	int		nboltz;		// number of Boltzmanns added
	int		i;			// loop counter

	// calculate number of Boltzmanns added together
	// check number of parameters
	// if error, return zero for f(p, x) & df/dp

	nboltz = (npar - 2) / 3;
	if (npar != (3 * nboltz + 2)) {
		for (i = 0; i < npar; i++) {
			dfdp[i] = 0.0;
		}
		return(0.0);
	}

	// initialize f(x) & df/dp

	y = p[0];
	dfdp[0] = 1.0;
	dfdp[1] = 0.0;
	drvforce = x - p[1];

	// loop through Boltzmann sums
	// add sums to f(x) & df/dp
	// calculate relevant df/dp at each sum

	for (i = 1; i <= nboltz; i++) {
		expfac = exp(-p[3 * i] * (x - p[3 * i + 1]));
		expprod = 1.0 / (1.0 + expfac);
		dfdp[1] += p[3 * i - 1] * expprod;
		dfdp[3 * i - 1] = drvforce * expprod;
		expprod *= expprod * expfac;
		dfdp[3 * i] = drvforce * p[3 * i - 1] * (x - p[3 * i + 1]) * expprod;
		dfdp[3 * i + 1] = -drvforce * p[3 * i - 1] * p[3 * i] * expprod;
	}
	y += drvforce * dfdp[1];
	return(y);
}
